#rmarkdown::render(“ClusteringCoseq.Rmd”)
On cherche à construire un réseau de régulation, sur les gènes qui ont été détectés comme différentiellement exprimés lors du traitement au KNO3 des plantes WT.
Etude de la quantification
data <- read.csv("quantifFiles/quantif_At.csv", h = T, sep = ",")
rownames(data) <- data$Gene
genes = which(!(grepl("__", rownames(data))))
not_quant = data[which((grepl("__", rownames(data)))), ]
data = data[genes, ]
Quant_genes = sum()
res = data.frame(No_features = as.numeric(as.vector(not_quant["__no_feature", grepl("R",
colnames(not_quant))])), Ambiguous = as.numeric(as.vector(not_quant["__ambiguous",
grepl("R", colnames(not_quant))])), Alignment_not_unique = as.numeric(as.vector(not_quant["__alignment_not_unique",
grepl("R", colnames(not_quant))])), QuantGenes = colSums(data[, grepl("R", colnames(data))]))
library(reshape2)
resMetlted <- melt(res)
p <- ggplot(resMetlted, aes(x = variable, y = value, fill = variable)) + geom_boxplot(alpha = 0.4) +
ggtitle("Nombre de reads par catégorie après quantification") + geom_jitter(width = 0.1,
alpha = 0.2)
p
Il semble qu’il y ait beaucoup de gènes ambigus pour htseq-count : cela se produit quant une coordonnée génomique chevauche plusieurs gènes. Peut être qu’on a de grands reads, et qu’ils chevauchent souvent les UTR des autres?
On refait l’analyse avec un autre pipeline (Cutadapt, Hisat2, Stringtie) et on observe la somme des counts par échantillons.
dataS <- read.csv("quantifFiles/quantifStringtie.csv", h = T, sep = ",")
rownames(dataS) <- dataS$gene_id
res$QuantGenesStringtie = colSums(dataS[, grepl("R", colnames(dataS))])
library(reshape2)
res <- melt(res)
p <- ggplot(res, aes(x = variable, y = value, fill = variable)) + geom_boxplot(alpha = 0.4) +
ggtitle("Nombre de reads par catégorie après quantification") + geom_jitter(width = 0.1,
alpha = 0.2)
p
On en conclue que Stringtie ne s’embête pas et compte les reads chevauchants pour chaque partie chevauchée, expliquant ce plus grand nombre de reads quantifiés.
On affiche la même chose pour mon htseq-count sur l’alignement de Hisat2 :
data <- read.csv("quantifFiles/quantifHtseqHisat.csv", h = T, sep = ",")
rownames(data) <- data$Gene
genes = which(!(grepl("__", rownames(data))))
not_quant = data[which((grepl("__", rownames(data)))), ]
data = data[genes, ]
Quant_genes = sum()
res = data.frame(No_features = as.numeric(as.vector(not_quant["__no_feature", grepl("R",
colnames(not_quant))])), Ambiguous = as.numeric(as.vector(not_quant["__ambiguous",
grepl("R", colnames(not_quant))])), Alignment_not_unique = as.numeric(as.vector(not_quant["__alignment_not_unique",
grepl("R", colnames(not_quant))])), QuantGenes = colSums(data[, grepl("R", colnames(data))]))
library(reshape2)
resMetlted <- melt(res)
p <- ggplot(resMetlted, aes(x = variable, y = value, fill = variable)) + geom_boxplot(alpha = 0.4) +
ggtitle("Nombre de reads par catégorie après quantification") + geom_jitter(width = 0.1,
alpha = 0.2)
p
On compare htseq-count avec l’option nonunique all, qui doit faire comme stringtie.
On essaie de faire un mapping plus stringent avec STAR en autorisant que deux mismatchs, mais on a des résultats très très similaires :
data <- read.csv("quantifFiles/quantif_At_Stringent.csv", h = T, sep = ",")
rownames(data) <- data$Gene
genes = which(!(grepl("__", rownames(data))))
not_quant = data[which((grepl("__", rownames(data)))), ]
data = data[genes, ]
Quant_genes = sum()
res = data.frame(No_features = as.numeric(as.vector(not_quant["__no_feature", grepl("R",
colnames(not_quant))])), Ambiguous = as.numeric(as.vector(not_quant["__ambiguous",
grepl("R", colnames(not_quant))])), Alignment_not_unique = as.numeric(as.vector(not_quant["__alignment_not_unique",
grepl("R", colnames(not_quant))])), QuantGenes = colSums(data[, grepl("R", colnames(data))]))
library(reshape2)
resMetlted <- melt(res)
p <- ggplot(resMetlted, aes(x = variable, y = value, fill = variable)) + geom_boxplot(alpha = 0.4) +
ggtitle("Nombre de reads par catégorie après quantification") + geom_jitter(width = 0.1,
alpha = 0.2)
p
La seule différence est qu’on a supprimé les outliers supérieurs dans la quatégorie des gènes quantifiés.
La bonne quantification
On a un soucis dans l’annotation : on a plusieurs features par transcrit! Les versions .1, .2, se chevauchaient tout le temps, problème repéré par visu dans igv. On choisit donc un seul transcrit par gène, et on re-quantifie.
data <- read.csv("quantifFiles/quantifNewAnnot.csv", h = T, sep = ",")
rownames(data) <- data$Gene
genes = which(!(grepl("__", rownames(data))))
not_quant = data[which((grepl("__", rownames(data)))), ]
data = data[genes, ]
Quant_genes = sum()
res = data.frame(No_features = as.numeric(as.vector(not_quant["__no_feature", grepl("R",
colnames(not_quant))])), Ambiguous = as.numeric(as.vector(not_quant["__ambiguous",
grepl("R", colnames(not_quant))])), Alignment_not_unique = as.numeric(as.vector(not_quant["__alignment_not_unique",
grepl("R", colnames(not_quant))])), QuantGenes = colSums(data[, grepl("R", colnames(data))]))
library(reshape2)
resMetlted <- melt(res)
p <- ggplot(resMetlted, aes(x = variable, y = value, fill = variable)) + geom_boxplot(alpha = 0.4) +
ggtitle("Nombre de reads par catégorie après quantification") + geom_jitter(width = 0.1,
alpha = 0.2)
p
data <- read.csv("quantifFiles/quantifGenes.csv", h = T, sep = ",")
rownames(data) <- data$Gene
genes = which(!(grepl("__", rownames(data))))
not_quant = data[which((grepl("__", rownames(data)))), ]
data = data[genes, ]
Quant_genes = sum()
res = data.frame(No_features = as.numeric(as.vector(not_quant["__no_feature", grepl("R",
colnames(not_quant))])), Ambiguous = as.numeric(as.vector(not_quant["__ambiguous",
grepl("R", colnames(not_quant))])), Alignment_not_unique = as.numeric(as.vector(not_quant["__alignment_not_unique",
grepl("R", colnames(not_quant))])), QuantGenes = colSums(data[, grepl("R", colnames(data))]))
library(reshape2)
resMetlted <- melt(res)
p <- ggplot(resMetlted, aes(x = variable, y = value, fill = variable)) + geom_boxplot(alpha = 0.4) +
ggtitle("Nombre de reads par catégorie après quantification") + geom_jitter(width = 0.1,
alpha = 0.2)
p
Homogénéité des réplicats : visualisation
annot <- read.csv("Code_for_RNAseq_CO2_N_Fr.csv", h = T, sep = ";")
conditions <- as.vector(unique(annot$Sample))
annot$condition <- substr(conditions, 1, nchar(conditions) - 1)
cond <- unique(substr(conditions, 1, nchar(conditions) - 1))
cond <- cond[grepl("At", cond)]
cor_diff <- c()
cor_same <- c()
scatter <- function(df, x, y) {
sp <- ggplot(df, aes(log(df[, x]), log(df[, y]))) + geom_bin2d(bins = 70) + scale_fill_gradient2() +
theme(legend.position = "none") + labs(x = x, y = y)
return(sp)
}
for (co in cond) {
print(co)
samples <- paste0("R", annot[annot$condition == co, ]$Code)
df <- data[, samples]
a <- scatter(df, samples[1], samples[2])
b <- scatter(df, samples[1], samples[3])
c <- scatter(df, samples[2], samples[3])
grid.arrange(a, b, c, ncol = 3)
cor_same <- c(cor_same, cor(df[, samples[1]], df[, samples[2]]), cor(df[, samples[1]],
df[, samples[3]]), cor(df[, samples[2]], df[, samples[3]]))
other <- paste0("R", annot[annot$condition != co & !grepl("Sl", annot$Sample),
]$Code)
for (o in other) {
cor_diff <- c(cor_diff, cor(data[, o], df[, samples[1]]), cor(data[, o],
df[, samples[2]]))
}
}[1] "At_AmbientCO2_HighNitrate_Fe"
[1] "At_AmbientCO2_LowNitrate_Fe"
[1] "At_ElevatedCO2_HighNitrate_Fe"
[1] "At_ElevatedCO2_LowNitrate_Fe"
[1] "At_AmbientCO2_HighNitrate_FeStarvation"
[1] "At_AmbientCO2_LowNitrate_FeStarvation"
[1] "At_ElevatedCO2_HighNitrate_FeStarvation"
[1] "At_ElevatedCO2_LowNitrate_FeStarvation"
Homogénéité des réplicats : significativité
On cherche maintenant à voir si ces plots montrent suffisamment corrélés pour montrer que les réplicats au sein d’une même condition sont suffisamment similaires. On peut déjà vérifier en plottant deux réplicats de conditions différentes : on voit que le scatter plot est plus dispersé que les réplicats biologiques.
cors = data.frame(`Different conditions` = cor_diff, `Same conditions` = cor_same)
library(reshape2)
cors <- melt(cors)
p <- ggplot(cors, aes(x = value, fill = variable)) + geom_density(alpha = 0.4) +
ggtitle("Distribution des corrélations entre les réplicats, biologiques ou non")
p
Nous pouvons donc en conclure que nous réplcats on bien été labelisés, et qu’ils sont plus homogènes au sein d’une même condition qu’entre différentes conditions de la combinatoire. Il semblerait que nous puissions passer à la suite : analyse stat pour l’expression différentielle.